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We systematically compare an event-by-event heavy-ion collision model to data from the Large 
Hadron Collider. Using a general Bayesian method, we probe multiple model parameters including 
fundamental quark-gluon plasma properties such as the specific shear viscosity 77 /s, calibrate the 
model to optimally reproduce experimental data, and extract quantitative constraints for all pa¬ 
rameters simultaneously. The method is universal and easily extensible to other data and collision 
models. 


I. INTRODUCTION 

Relativistic heavy-ion collisions produce a hot, dense 
phase of strongly-interacting matter commonly known 
as the quark-gluon plasma (QGP), which rapidly ex¬ 
pands and freezes into hadrons jl 7 . Since the QGP 
is not directly observable—only final-state hadrons are 
detected—present research seeks to quantify the funda¬ 
mental properties of the QGP, such as its transport coef¬ 
ficients and the nature of the initial state, through com¬ 
parisons of experimental measurements to computational 
model calculations. 

Computational models must take a set of input pa¬ 
rameters including the physical properties of interest, 
simulate the full time-evolution of heavy-ion collisions, 
and produce outputs analogous to experimental measure¬ 
ments. The true values of the physical properties are ex¬ 
tracted by calibrating the input parameters so that the 
model output optimally reproduces experimental data. 
This generic recipe is called “model-to-data comparison”. 

Notably, the QGP shear viscosity to entropy density 
ratio r]/s has been constrained by comparing anisotropic 
flow coefficients v n between model and experiment. Ex¬ 
plicit calculation of r]/s directly from QCD is not yet 
feasible, and while there is a conjectured lower bound 
77/s > 1/47T ~ 0.08 from AdS/CFT holography [8], 
model-to-data comparison is the most attractive option 
for determining the optimal input parameter value and 
corresponding uncertainty. To this end, previous studies 
used viscous relativistic fluid dynamics and hybrid trans¬ 
port models to compute v n at several values of 77/s, then 
chose the value which most closely matched experimental 
v n . A variety of complementary calculations have con¬ 
strained 77/s to an approximate range of 0 . 08 - 0.20 9 12]. 

However, 77/s is not the only model input parame¬ 
ter: many other parameters remain poorly determined, 
e.g. the hydrodynamic thermalization time To and initial 
conditions; and models often have non-physical nuisance 
parameters that nonetheless should be tuned to optimal 
values. The flow coefficients v n are only a small subset 


of all QGP observables: models must also describe basic 
quantities such as the charged-particle multiplicity and 
transverse-momentum distributions. 

Recent work [13 moved toward a more global anal¬ 
ysis of multiple model parameters and observables, but 
encountered practical limitations attempting to simulta¬ 
neously tune these free parameters. In general, input 
parameters correlate among each other and contribute 
to multiple observables, so they cannot be constrained 
independently. 

Algorithms such as Markov chain Monte Carlo 
(MCMC) can rigorously explore this type of complex 
high-dimensional parameter space, but require a very 
large number of model evaluations—often thousands or 
millions, depending on the problem at hand. Heavy- 
ion collision models may run for several hours, so a di¬ 
rect MCMC approach is intractable. The situation is 
exacerbated when studying event-by-event fluctuations 
as opposed to average quantities: while event-averaged 
models save computation time by using a smooth initial 
condition and single hydrodynamic calculation, event-by¬ 
event models have realistic, fluctuated initial conditions, 
each of which requires its own hydrodynamic treatment. 
Many thousands of complete events are necessary at each 
point in parameter space to capture event-by-event fluc¬ 
tuations. 

These limitations may be overcome through a modern 
Bayesian method for analyzing computationally expen¬ 
sive models mm- a set of salient model parameters is 
chosen for calibration—the set should include any funda¬ 
mental physical properties of interest—and the model is 
evaluated at a relatively small 0(1O 2 ) number of points. 
Those points are then interpolated with a Gaussian pro¬ 
cess emulator El to provide a continuous picture of 
the parameter space. The emulator acts as a fast sur¬ 
rogate to the full model: it predicts model output at 
arbitrary points in parameter space with negligible com¬ 
putational cost. This effectively removes most practical 
barriers and enables parameter calibration through stan¬ 
dard techniques such as MCMC. 
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Emulators have been successfully used to study a wide 
range of physical systems, including galaxy formation 
[18; and heavy-ion collisions HMD- Reference m cali¬ 
brated a hydrodynamic model to identified particle spec¬ 
tra from the Relativistic Heavy Ion Collider (RHIC) and 
extracted constraints on rj/s and several initial state pa¬ 
rameters. However, this study used an event-averaged 
initial condition model, limiting its ability to investigate 
event-by-event fluctuations. 

In this work, we apply Bayesian methodology to a full 
event-by-event heavy-ion collision model. We calibrate to 
multiplicity and flow data from the Large Hadron Col¬ 
lider (LHC) and constrain the shear viscosity r]/s along 
with other hydrodynamic and initial condition parame¬ 
ters. The analysis framework handles arbitrary numbers 
of inputs and outputs, systematically calculates quantita¬ 
tive constraints on all inputs simultaneously, and quickly 
evaluates the efficacy of physical models. 

II. MODEL 

State-of-the-art heavy-ion collision models simulate 
QGP spacetime evolution in several stages HQumcEa- 

m 

1. an initial condition model describes the initial state 
and non-equilibrium dynamics until QGP forma¬ 
tion, 

2. viscous relativistic hydrodynamics calculates the 
dynamical expansion of the hot and dense QGP 
medium including the phase transition to a hadron 
gas, 

3. then a particlization model converts the system into 
a microscopic ensemble of hadrons, 

4. and finally a Boltzmann transport model calculates 
hadronic rescattering and decays. 

In this work, we opt for a mature, well-tested set of event- 
by-event models m with an established track record of 
describing diverse RHIC and LHC data jT(). [28. [29 . This 
choice will permit direct comparison between existing re¬ 
sults and the outcome of the following systematic model- 
to-data comparison. We emphasize, however, that the 
methodology in this paper can easily be applied to any 
set of models and corresponding data. 

A. Initial conditions 

Initial condition models provide the outcome of the col¬ 
lision’s pre-equilibrium evolution at the hydrodynamic 
thermalization time, approximately 0.5 fm/c. Some 
models explicitly calculate pre-equilibrium dynamics m 
starting from the initial state of the collision; others skip 
this time frame and generate initial conditions directly 
at the thermalization time [3TH33] . 


We select two of the most widely used models in the 
latter category: the Monte Carlo Glauber [32] and Monte 
Carlo KLN f3L models. Although more sophisticated 
models were recently introduced [30j [33] , both Glauber 
and KLN provide reasonable event-by-event initial con¬ 
ditions with well-understood behavior and a broad basis 
of published results. 

B. Hydrodynamics 

The initial condition furnishes the hydrodynamic 
stress-energy tensor at the thermalization time 

To- Viscous hydrodynamics then solves the conservation 
equations 

d^ v = 0 ( 1 ) 

where 

= (e + P)u»u v - Pg» v + tt^; (2) 

e, P, and u M are the energy density, pressure, and flow 
velocity of the fluid; g^ v is the metric tensor; and tt^ is 
the shear stress tensor. An equation of state 

P = P(e) (3) 

closes the system of hydrodynamic equations and is usu¬ 
ally provided by a parametrization of lattice QCD calcu¬ 
lations. 

We employ an improved version of VISH2+1 [34], a 
stable, extensively tested implementation of boost-invariant 
viscous hydrodynamics that was recently updated to handle 
fluctuating event-by-event initial conditions [27]. VISH2+1 
uses the prevalent s95 partial chemical equilibrium equation 
of state [55] . 

C. Particlization 

As the hydrodynamic medium expands and cools, it un¬ 
dergoes a transition from a deconfined QGP to a hot and 
dense hadronic system. At this point it’s advantageous to 
switch to a microscopic transport model, for such models 
naturally account for the system’s increasing viscosity, non¬ 
equilibrium break-up, and eventual freeze-out. Thus, a parti¬ 
clization model converts the fluid into a microscopic ensemble 
of hadrons once the fluid cools to a pre-specified switching 
temperature, typically just below the QCD transition tem¬ 
perature T c ~ 165 MeV. The model generates particles by 
sampling the Cooper-Frye formula [36] 

E ^P = f d3(7 V’ ( 4 ) 

where fi is the distribution function for particle species i, 
is the four-momentum, and the integral is taken over the 
isothermal spacetime hypersurface a defined by the switching 
temperature. 

We use a recent hypersurface sampler designed to couple 
with VISH2+1 [27l 137] , 
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D. Hadronic phase 

After particlization, the medium continues to expand as 
an interacting hadron gas (e.g. scatterings and decays). A 
hadronic “afterburner” calculates these interactions through 
the Boltzmann equation 

^=&(z,p), (5) 

where fi is the distribution function and Ci is the collision 
kernel which contains all possible hadronic interactions for 
particle species i. Particles emerging from the afterburner 
are analogous to particles streaming into an experimental de¬ 
tector. 

We adopt Ultra-relativistic Quantum Molecular Dynamics 
(UrQMD) [38, [39] as an afterburner. 


E. Postprocessing 

The full event-by-event model is executed 0( 10 4 ) times for 
each set of input parameters, yielding 0( 10 7 ) events in to¬ 
tal. Events are binned into centrality intervals and the raw 
event data are postprocessed into physical observables for di¬ 
rect comparison with experiment. In this analysis we calcu¬ 
late the centrality dependence of several standard observables: 
the average charged-particle multiplicity (iV c h) and multi¬ 
particle flow cumulants v n {2k}. Note that the method triv¬ 
ially extends to arbitrary numbers and types of observables— 
all that’s required is a model calculation and a corresponding 
experimental measurement. 

Flow cumulants v n {2k} are defined as the 2/c-particle cor¬ 
relation function of the nth-order azimuthal anisotropy. For 
example, the two-particle cumulant is 

v n { 2} 2 = (6) 

where <fii is the azimuthal angle of the transverse momentum 
of particle i and the average is over all distinct pairs of par¬ 
ticles i,j. The two-particle cumulant is also approximately 
equal to the root-mean-square of the full v n distribution m- 

v n {2) ~ v / (wU- (7) 

We compute two-particle cumulants for elliptic and triangu¬ 
lar flow ^ 2 {2}, ^ 3 { 2 } using the direct Q-cumulant method 
ED- Higher-order cumulants are currently out of reach due 
to insufficient quantities of events. 

Postprocessed observables are compared to corresponding 
experimental results recently measured by the ALICE exper¬ 
iment at the LHC for Pb-Pb collisions at ^/snn = 2.76 TeV 
M • All observables are subjected to the same kinematic cuts 
as the ALICE detector, namely charged particles with \rj\ < 1 
and 0.2 < pr < 3.0 GeV. 


III. EMULATOR 

This section constructs a Gaussian process (GP) emulator 
to act as a surrogate for the full event-by-event model. The 
strategy is to evaluate the model on a carefully chosen set of 
input parameter points, then use a GP to interpolate the pa¬ 
rameter space. Unlike alternative interpolation schemes such 


as splines or polynomial interpolation, a GP emulator pro¬ 
vides a probability distribution at each point in parameter 
space, hence, it not only predicts the output of the model 
at arbitrary points in parameter space, but also quantifies 
the uncertainty of its predictions. Further, GPs are non- 
parametric interpolators, i.e. they do not require an assumed 
functional form for the underlying model. These features are 
essential for emulation of computer codes. 


A. Gaussian processes 


This subsection summarizes the theory of Gaussian process 
emulators as detailed in Chap. 2 of HU. 

A Gaussian process (GP) is defined as a collection of ran¬ 
dom variables, any finite number of which have a joint Gaus¬ 
sian distribution. A GP may be thought of as a stochastic 
function /(x) which maps n-dimensional input vectors x to 
normally distributed outputs y. It is fully specified by a mean 
function /i(x) which gives the mean of / at input point x and 
a covariance function cr(x, xQ which provides the covariance 
of / between a pair of points x, xb 

As a concrete example, let xi be an input point and 
yi = /(xi) be the output of the GP at xi; then y\ has a 
normal distribution with mean /i(xi) and variance cr(xi,xi): 

3/1 ~ A/'(p(xi),cr(xi,Xi)). (8) 


Now if x 2 is another input and 7/2 = /(x 2 ) is the correspond¬ 
ing output, yi and y 2 have a bivariate normal distribution 


yi i-A/UTGi CU’U T- X2 b 

y 2 J |_ W x 2)/ ’ V7 ( x 2, Xi ) Cr(x 2 ,X 2 )y 


(9) 


In general, the set of m random output variables 
y — { 2 / 1 ,... , 2 / m } =1 f(X) corresponding to input points 
X = {xi,...,x m } have a multivariate normal distribution 

y (10) 

where 

M = K X ) = M X 0> M(x 2 ), • • • , /t(Xm)} (11) 

is the m -dimensional mean vector from applying the mean 
function to each input, and 

( cr(xi,xi) ••• <r(xi,x m )\ 

: : ( 12 ) 

<j(x m ,xi) ••• a (x m ,x m )/ 


is the m x m covariance matrix from applying the covariance 
function to each pair of inputs. 

In practice, the mean function is often set to zero, since 
the mean of a distribution can always be subtracted off. The 
covariance function must be carefully chosen, for it controls 
the degree of similarity between pairs of points. A standard 
choice is the squared-exponential function 

g(x,x') =exp( L X 2 ^ C I ), (13) 


where t is a characteristic length scale. Notice that nearby 
points are strongly correlated (awl) while distant points 
approach independence (a —> 0). This implies that the GP is 
smooth, i.e. nearby input points produce similar outputs. 
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Just as we can sample random numbers from a distribution, 
we can draw random functions from a GP. We choose a set of 
test points X* (the reason for the subscript * will become clear 
in a moment), calculate the covariance matrix E as cr(X*, X*), 
and generate multivariate normal samples from A/"(0, E). We 
can then plot the input-output points as smooth curves, as in 
the top panel of Fig. [l] 

Of course, simply generating random functions is not par¬ 
ticularly useful—we want to use a GP to interpolate a com¬ 
puter model. Suppose we have a model which takes a vector 
of input parameters x and produces an output y according 
to some unknown GP /(x); for example, / could be a hydro- 
dynamic model with input parameters x = (to, 77 /s) and the 
output could be elliptic flow V 2 - We choose a set of train¬ 
ing points X , run the model at each point, and observe a set 
of outputs y. Now, instead of completely random functions, 
we desire functions which pass through (interpolate) all the 
training points (X, y). This is achieved by conditioning the 
GP on the training data to yield a predictive distribution for 
y at any input point x. Recalling the test points X*, the pre¬ 
dictive distribution for the corresponding outputs y* is the 
multivariate normal distribution 

y* ~ V(m,s), 

fi = a(X*,X)a(X,X )~ 1 y, (14) 

S = cr(X*,X*) - £r(X,,X)<T(X,X) _1 0 -(X,X,). 

See the bottom panel of Fig. [TJfor an example of conditioning 
a GP on one-dimensional training data. 

We emphasize that the prediction y* is not constant, but 
a probability distribution for the model outputs at X*. As 
demonstrated in Fig. [l] the predictive distribution is narrow 
when near the training points and wide when far away, hence, 
it reflects the true state of knowledge of the interpolation. 
This is accomplished without assuming a parametric form for 
the model—we must only assume that the model is a GP with 
a specified covariance function. 
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FIG. 1. Top: Random functions drawn from a Gaussian 
process using a squared-exponential covariance function with 
length scale t = 1. Bottom: Functions drawn from a GP 
conditioned on the training points indicated by dots. In both 
plots, the dashed line represents the GP mean and the grey 
band is twice the GP standard deviation (roughly 95% confi¬ 
dence interval). 


B. Computer experiment design 

The full event-by-event model is to be evaluated on a set of 
m training points X & {xi,..., x m }, where each x* is an n- 
dimensional vector of input parameters, so X may be viewed 
as an m x n design matrix. This subsection details the choice 
of input parameters and their distribution in parameter space. 

For the present study, we choose a set of n — 5 input pa¬ 
rameters 

x= (Norm, I.C. param, to, 77 /s, r^) (15) 

where 

• Norm is the overall normalization factor, a multiplica¬ 
tive constant that determines how much entropy is de¬ 
posited in the initial condition. 

• I.C. param is a parameter specific to each initial condi¬ 
tion model. For the Glauber model the parameter is a , 
which controls how entropy is distributed to wounded 
nucleons and binary collisions; for the KLN model it 
is A, a dimensionless exponent in the saturation scale 
parametrization. Both are related to the centrality de¬ 
pendence of multiplicity. 


• To is the QGP thermalization time and the starting time 
for hydrodynamic evolution. 

• 77 /s is the shear viscosity to entropy density ratio of 
the QGP, assumed to be fixed throughout the hydro- 
dynamic evolution stage. 

• Ttt is the shear stress relaxation time, which dictates 
how quickly the hydro medium relaxes to the Navier- 
Stokes limit. Since the relaxation time is a function of 
the shear viscosity and temperature and thus cannot be 
tuned explicitly, we use the coefficient k n in the relation 
t^ =£ 6k n rj /(sT) as a tunable parameter. 

We set intentionally large ranges for each parameter, summa¬ 
rized in Table [I] In this work, we fix several auxiliary parame¬ 
ters to reasonable defaults: nucleons are assumed to be disks 
with size determined by the inelastic nucleon-nucleon cross 
section <tnn, and the hydro to micro switching temperature 
is set to 165 MeV, just below the equation of state transi¬ 
tion temperature. However, the method can handle arbitrary 
numbers of parameters provided the sample size is sufficiently 
large. 

The training points X = {xi,... ,x m } must be chosen to 
simultaneously optimize emulator accuracy and computation 
time. Perhaps the most obvious design strategy is a uni¬ 
form grid (factorial design), e.g. k evenly-spaced points in 
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TABLE I. Input parameter ranges for the Glauber (Gib) 
and KLN initial condition models and for the hydrodynamic 
model. 


Parameter 

Description 

Range 

Gib Norm 

Overall normalization 

20-60 

Gib a 

Wounded nucleon / binary coll. 

0.05-0.30 

KLN Norm 

Overall normalization 

5-15 

KLN A 

Saturation scale exponent 

0.1-0.3 

TO 

Thermalization time 

0.2-1.0 fm 

p/s 

Specific shear viscosity 

0-0.3 

k^ 

Shear relaxation time coeff. 

0.2-1.1 


each dimension. Unfortunately, this leads to a total sample 
size m — k n which even for a modest k = 10 and n — 5 is 
intractably large. 

A popular algorithm for generating efficient design points 
is maximin Latin hypercube sampling ^S] • This method pro¬ 
duces space-filling randomized designs with several desirable 
properties: 

• The minimum distance between points is maximized , 
thus avoiding large gaps and tight clusters. 

• Projections of the design into lower dimensions are uni¬ 
formly distributed. 

Figure [2] illustrates these traits. A Latin hypercube design 
with a relatively small sample size provides an efficient scaf¬ 
folding of parameter space for interpolation by a GP emulator. 
As a general rule of thumb, a sample size m ~ lOn yields ac¬ 
ceptable interpolation accuracy m and is a common choice 
for an initial experiment with limited computation time, how¬ 
ever there is no harm in a larger sample. 

We use a 256 point Latin hypercube design across the n = 5 
input parameters; Fig.[2]shows a two-dimensional projection. 
At each design point, we have executed the event-by-event 
model 0( 10 3 ) times in each of six centrality bins 0-5%, 10- 
15%, ..., 50-55%, for both the Glauber and KLN models, 
yielding G( 10 7 ) events in total. Two design points that were 
very near to the edge of the design space gave non-physical 
results and have been discarded, so the operational design has 
m — 254 points. 

Figure |3] shows the postprocessed observables (iV c h), M2}, 
V 3 {2} as a function of centrality for each point in the design. 
The results have a broad distribution which is a direct result of 
the wide ranges of input parameters. There is some statistical 
error present in v 3 {2} due to insufficient quantities of events. 

Note that these results constitute the training data for the 
GP emulator, not any kind of best-fit. 


C. Multivariate output 

Gaussian processes are fundamentally scalar functions, but 
computer models often produce multivariate output. In gen¬ 
eral, the model takes the mxn design matrix X and computes 
an mxp output matrix Y. The present event-by-event model 
has p = 18 outputs (three observables each in six centrality 
bins). 
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FIG. 2. The Latin hypercube experiment design projected 
into the ( 77 /s, to) dimensions. All other parameters also vary 
across the design, so the points that appear very close in 
the projection are not necessarily close in the full-dimensional 
space. The edge histograms show the distributions flattened 
into one dimension; note that they are space-filling and ap¬ 
proximately uniform. 

An obvious workaround is to use independent GP emulators 
for each of the p outputs, however, this would neglect corre¬ 
lations and quickly become unwieldy for higher-dimensional 
output. Instead, we decompose the outputs into orthogonal 
linear combinations called principal components (PCs) and 
emulate each transformed PC. The PCs are uncorrelated by 
construction and can also be used to reduce the dimension¬ 
ality of the output space. Figure shows an example PC 
decomposition. 

To calculate the PCs, we first subtract the mean of the 
output data Y so that each column has mean zero, then com¬ 
pute the eigendecomposition of the sample covariance matrix 
Y T Y: 

Y J Y = UAU\ (16) 

where U is an orthogonal p x p matrix containing the eigen¬ 
vectors of Y T Y and A is diagonal containing the eigenvalues 
Ai,...,A p in non-increasing order. U now defines a linear 
transformation which “rotates” the output data Y into PC 
space: 

Z = y/rnYU, (17) 

where Z is an m x p matrix (same shape as 7) of the trans¬ 
formed PCs. The eigenvalues A* represent the variance ex¬ 
plained by principal component i; since they are sorted in 
non-increasing order, the fraction of the variance explained 
by the first q < p PCs is 

V(q) = HA/ (18) 

2^=i 

Often, the first few PCs describe most of the variance, as 
demonstrated for the present data in Fig. [5] Hence we can 
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FIG. 3. Model calculations from Glauber (top, blue) and KLN (bottom, green) initial conditions. Each plot has 254 lines 
corresponding to the 254 Latin hypercube design points. From left to right: average charged-particle multiplicity (iV c h), elliptic 
flow two-particle cumulant ^ 2 {2}, and triangular flow two-particle cumulant u 3 {2}. Data points are experimental measurements 
from ALICE [42]. 


construct a reduced-dimension transformation with minimal 
loss of precision by choosing q < p so that V(q) satisfies 





FIG. 4. Principal component decomposition of the observ¬ 
ables yj (iV c h), ^ 2 ( 2 } for the Glauber model in 20-25% cen¬ 
trality. Each data point represents a model calculation and 
the edge histograms show the approximate normal distribu¬ 
tion of each observable. Arrows represent the PC vectors with 
lengths proportional to the explained variance. 


some threshold (e.g. V(q) > 0.99) and taking only the first q 
columns of U : 


Zq — yJrnYUq , 


(19) 


where Z q is now an m X q matrix. 

We may now use q independent GP emulators for each of 
the columns of Z q . GPs are conditioned on the design X 
according to Eq. (14) and predict the PCs Z* at arbitrary 
test points A* which are then transformed back to physical 
space as 


F, = 4z,U t . (20) 

\/m 

There is an important caveat for principal components: the 
original data Y must have a multivariate normal distribu¬ 
tion for the transformed PCs Z to be uncorrelated. There is 
no guarantee that a particular model will produce normally- 
distributed outputs so this must be verified on a case-by-case 
basis. For the present event-by-event model we perform the 
following steps: 


1. Assess the normality of each observable (iV c h), ^ 2 ( 2 }, 
u 3 {2}. While the flow cumulants are approximately 
normal without modification, we take the square root 
of multiplicity yj{N c h) to obtain a normal distribution, 
as shown in Fig. [4] 

2. Divide each observable by its corresponding experimen¬ 
tal value from ALICE !42] . This converts everything to 
unit less quantities of order one. 

3. Multiply each observable by a manually specified 
weight factor, ratios 1.2 : 1.0 : 0.6 for observables 
(iVch) • ^ 2 ( 2 } : u 3 {2}. These subjective weights encour¬ 
age the model to fit more strongly to more fundamental 
observables, e.g. we prefer a model that describes N c h 
and V 2 at the expense of u 3 rather than fitting u 3 with 
incorrect N c h. The weights will be discussed further in 
the results, Sec. IE! 
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FIG. 5. Fraction of the variance V(q) explained by the first 
q principal components for Glauber (blue circles) and KLN 
(green triangles), q — 5 explains approximately 99% of the 
total variance, a significant reduction from the original 18 
dimensions. 


4. Concatenate the unitless, weighted data into a 254 x 18 
matrix Y, where each row corresponds to a design point 
and each column to an observable and centrality bin. 

5. Subtract the mean of each column and transform Y 
into principal components Z q with q — 5 PCs, retaining 
over 99% of the variance as shown in Fig. [5] The PC 
transformation matrices U q , shown in Fig. pp reflect the 
natural correlations among observables, for example all 
observables are correlated in the first PC (~75% of total 
variance), while iV c h is anti-correlated with the v n in the 
second PC (~20% of total variance). 

We invert these steps to transform PCs back to physical space. 

In practice, the covariance method for computing principal 
components is prone to numerical error, so a more robust 
algorithm using the singular value decomposition (SVD) is 
preferred. The SVD of the data Y is 

Y m VDW J (21) 


where V, W are orthogonal matrices containing the so-called 
left- and right-singular vectors of Y and D is diag onal con¬ 
taining the singular values. Inserting (21) into (16) yields 


Y J Y = WD W J = UXU\ 


( 22 ) 


hence the singular values D are the square root of the eigen¬ 
values A and the right singular vectors W are the eigenvectors 


U. 


D. Constructing and validating the emulator 


We emulate the model by conditioning independent Gaus¬ 
sian processes on each of the principal components Z q and 
the input design X according to Eq. (14). Model outputs 


inevitably include statistical noise, i.e. we cannot compute 
y — /(x) exactly, only y — /(x) + e where e is Gaussian noise. 
This is accounted for by adding a noise term to the diagonal 
of the covariance matrix: 


r(x,x') - 5 - <t(x, x') + al5 x 


Glauber KLN 
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FIG. 6. Visualization of the principal component transfor¬ 
mation matrices U q for Glauber (left) and KLN (right). The 
numerical values of each matrix element are annotated and 
color-coded, where darker red indicates more positive values, 
darker blue indicates more negative, and grey indicates zero. 


where o^ is the variance of the noise and 5 XX / is a Kronecker 
delta. Effectively, the noise term relaxes the requirement that 
the GP must pass exactly through each training point. 

We use a squared-exponential covariance function with a 
noise term: 


/ /\ 2 
cr(x, x ) = o’ q P exp 


-E 

L k =1 


{xk - 4 ) 2 

K 


+ 


n "xx' i 


(23) 


where Oq P is the overall variance of the GP and ik is the 
characteristic length scale for dimension k. These hyperpa¬ 
rameters (crop, cr n , £k) are not known a priori and must be 
estimated from the training data, however, in the present case 
predictions appear to be relatively insensitive to the precise 
choice of hyperparameters, as will be demonstrated promptly. 
For details about the selection of hyperparameters see the |Ap-| 
|pendix| 

As with any interpolation scheme, the GP emulator must 
be validated to ensure it faithfully predicts model output. In 
other words, given an arbitrary test point x*, the GP pre¬ 
diction at x* should agree (within its uncertainty) with an 
explicit computation at x*. To this end, we have generated a 
separate 64-point Latin hypercube validation design X *, eval¬ 
uated the full event-by-event model at each validation point 
just as for the training design V, and predicted the model 
outputs at V* using the GP emulator. 

Figure [7] validates that the emulator does indeed faithfully 
predict the model. Recall that the emulator provides prob¬ 
ability distributions of finite width, so it need not predict 
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FIG. 7. Validation of the Gaussian process emulator for the Glauber model. Each plot shows emulator predictions against 
explicit calculations for the 64 validation design points in centrality bins 0-5% (green circles), 20-25% (orange triangles), and 
40-45% (purple squares). The x-value of each data point is the emulator prediction with 2cr (95%) horizontal error bars, the 
y -value is the explicit calculation with 2cr (95%) vertical error bars, and the diagonal grey line represents y = x. 


every validation point exactly—in fact, in the ideal case the 
residuals would have a normal distribution with mean zero. 
Most of the uncertainty visible in Fig. [7] is actually due to 
the statistical noise in the flow cumulants, especially ^ 3 { 2 }. 
The emulator accurately accounts for the noise present in the 
underlying data. 


IV. CALIBRATION 

With the validated Gaussian process emulator in hand, it 
may be used as a fast surrogate to the full event-by-event 
model for calibration. Calibration means tuning the model 
input parameters so that the output optimally agrees with 
experimental data and in the process extracting probability 
distributions for each parameter. Recall the input parameters 
are 

x= (Norm, I.C. param, To, 77 /s, k n ). 

Presumably there exists a true set of parameters x*; the task 
now is to find the probability distribution of x* given the 
training data (V, Y) and experimental measurements y exp . 
This distribution may be framed in terms of Bayes’ theorem 
as 

P(x*|V,y,y eX p) (x P(X, Y, y exp |x*)P(x*) (24) 

where 

• P(x*) is the prior probability which embodies initial 
knowledge of x*; 

• P(X, Y, y exp |x*) is the likelihood: the probability of 
observing (V, V, y exp ) given a proposed value of x*; and 

• P(x*|X, Y, y exp ) is the posterior probability for x* 
given the observations (X, Y, y exp ). This is the prob¬ 
ability distribution we wish to construct. 

In general Bayes’ theorem has a normalization constant which 
has been omitted since we are only concerned with relative 
probabilities. 

The remainder of this section applies the methodology from 
EHE] to calibrate the model and determine the posterior 
probability for x*. 


A. MCMC 

The workhorse of any Bayesian calibration is Markov chain 
Monte Carlo (MCMC), a powerful and flexible method for 
directly sampling the posterior probability. Perhaps the most 
common version is the Metropolis-Hastings algorithm, which 
generates a random walk through parameter space by accept¬ 
ing or rejecting steps based on the posterior probability. For a 
large number of steps the samples of the random walk equili¬ 
brate to the posterior distribution. We use the affine-invariant 
ensemble sampler for MCMC 55, H], an alternative algo¬ 
rithm that uses a large ensemble of interdependent walkers. 
Ensemble sampling notably has a much shorter autocorrela¬ 
tion time than Metropolis-Hastings sampling and hence con¬ 
verges more quickly to the equilibrium distribution. 

The MCMC algorithm samples proposal points x* and cal¬ 
culates the posterior probability at each point via Bayes’ the¬ 
orem. We place a non-informative flat prior on x*, that is, 
the prior probability is constant within the design range (Ta¬ 
ble [I]} and zero outside. We evaluate the likelihood in principal 
component space: 

-P(z e xp|x*) oc exp| - (z* z exp ) T F z (z* z eX p)^, (25) 

where z exp is the PC transform of the experimental data, 
z* is the emulator prediction of the PCs at x*, and 5% is 
the covariance (uncertainty) matrix on the PCs assuming 
normally-distributed errors. Given the flat prior, the pos¬ 
terior P(x* |z exp ) is equal to the likelihood within the design 
range and zero outside. 

There are a number of sources of uncertainty including 
experimental statistical and systematic error, model statis¬ 
tical and systematic error, and emulator uncertainty. In the 
present study, we do not attempt to precisely account for 
each contribution, for this would inevitably require dubious 
assumptions about systematic error correlations and the un¬ 
known error of the model. We assume a simple fractional 
error on the principal components, i.e. the covariance matrix 
is 


E z = diag(cr^ z eX p), 


(26) 
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FIG. 8. Posterior marginal and joint distributions of the calibration parameters for the Glauber model. On the diagonal are 
histograms of MCMC samples for the respective parameters, on the lower triangle are two-dimensional scatter histograms of 
MCMC samples showing the correlation between pairs of parameters, and on the upper triangle are approximate contours for 
68%, 95%, and 99% confidence regions along with a dot indicating the median. 


where <j 2 z is a manually set constant, cf z = 0.06 for the present 
study, to account for the typical experimental error of 3-5% 
[42] plus some additional uncertainty. While this is itself a 
rough assumption, it is perhaps no worse than the alternative, 
since experimental systematic errors are typically estimated 
percentages themselves and the principal component trans¬ 
formation automatically includes natural correlations among 
observables. The primary goal of this study is to develop and 


test a model-to-data comparison framework; details such as 
the precise treatment of uncertainties can be improved later. 

We run 0( 10 6 ) MCMC steps to allow the chain to equili¬ 
brate, discard these “burn-in” samples, then run 0( 10 7 ) steps 
to generate the posterior distribution. 
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FIG. 9. Same as Fig. [ 8 ] for the KLN model. 


B. Results 


The primary MCMC calibration results are presented in 
Figs. HI and [9] for the Glauber and KLN models, respectively. 
These are visualizations of the posterior probability distribu¬ 
tions of the true parameters x* , including the distribution of 
each individual parameter and all correlations. The diagonal 
histograms show the marginal distributions for each param¬ 
eter (all other parameters integrated out); the lower-triangle 
plots are two-dimensional scatter histograms of joint distribu¬ 


tions between pairs of parameters, where darker color denotes 
higher probability density; and the upper triangle has contour 
plots of the same joint distributions, where the contour lines 
enclose the 68 %, 95%, and 99% confidence regions. 

A wealth of information may be gained from these poste¬ 
rior visualizations; the following highlights some important 
features. 

Focusing on the Glauber results in Fig. [ 8 j we see the shear 
viscosity 77 /s (fourth diagonal plot) has a narrow approxi¬ 
mately normal distribution located near the commonly quoted 
value 0.08. As expected, rj/s is tightly constrained by exper- 
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FIG. 10. Random realizations of the calibrated posterior for Glauber (top, blue) and KLN (bottom, green) initial conditions. 
Similar to Fig. [3] except the lines are posterior emulator predictions instead of explicit prior calculations. 


imental flow data. Going across the fourth row, we observe 
nontrivial correlations among 77 /s and other parameters, for 
example, 77 /s and the hydrodynamic thermalization time To 
are negatively correlated (fourth row, third column). As To 
increases, the medium expands as a fluid for less time, so less 
flow develops, and viscosity must decrease to compensate. 

Both To and normalization (third and first diagonals) 
have broad distributions without strong peaks, and they are 
strongly-correlated (third row, first column). This is because 
the hydrodynamic model is boost-invariant and lacks any pre¬ 
equilibrium dynamics, so To is effectively an inverse normaliza¬ 
tion factor. The joint distribution shows a narrow acceptable 
band whose shape is governed by the inverse relationship. 

The wounded nucleon / binary collision parameter a (sec¬ 
ond diagonal) has a roughly-normal distribution located near 
the typical value 0.12. It is mainly related to the slope of mul¬ 
tiplicity vs. centrality and hence has a nontrivial correlation 
with normalization and To, e.g. we can decrease the normal¬ 
ization to the lower end of its distribution provided we also 
increase a to compensate. 

Meanwhile, the shear stress relaxation time coefficient k n 
(fifth diagonal) has an almost flat distribution and its joint 
distributions show no correlations. Evidently, this parameter 
does not influence flow coefficients or multiplicity. 

The KLN results in Fig. [9] generally exhibit wider, less nor¬ 
mal distributions than Glauber. This suggests that KLN is 
somewhat less flexible than Glauber, so its overall behavior 
is relatively insensitive to the specific values of input param¬ 
eters. 

The shear viscosity 77 /s has a narrow, irregular distribution 
covering the common value 0.20. As with Glauber, 77 /s has 
a negative correlation with to, there is a strong inverse rela¬ 
tionship between normalization and To, and k n has no effect. 
The KLN parameter A has a flat marginal distribution, but 
there are strongly excluded regions in the joint distributions 
with normalization and To. This appears to be the same effect 
as observed with Glauber a, except the dependence on A is 


significantly weaker. 

The posteriors may be validated by drawing samples from 
the calibrated distributions and visualizing the correspond¬ 
ing emulator predictions: if the model is correct and properly 
calibrated, the posterior samples will be close to experimental 
measurements. Figure [To| confirms—for the most part—that 
the posteriors are indeed tightly clustered around the data 
points. Visualizations such as this will always have some un¬ 
certainty since samples are drawn from the full posterior, how¬ 
ever, Fig. |10| has markedly narrower clusters than Fig. [3] in 
which the input parameters varied across their full ranges and 
were not tuned to match experiment. 

As shown in the top row of Fig. | 10 | the Glauber model 
nearly fits the centrality dependence of all the present observ- 


Glauber 



r//s 


FIG. 11. Comparison of posterior distributions of 77 /s for 
Glauber (blue) and KLN (green). These are the same his¬ 
tograms as in Figs. [8]and[9j expanded and placed on the same 
axis. The vertical grey lines indicate the common values 0.08 
for Glauber and 0.20 for KLN [25], 47]. 
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TABLE II. Quantitative summary of posterior distributions. For each parameter, the previous estimate EH El HE] , mean, 
median, and confidence intervals are given. 



Parameter 

Prev. est. 

Mean 

Median 

68 % C.I. 

95% C.I. 

99% C.I. 


Norm. 

57 

48.9 

49.0 

41.6-56.4 

36.5-59.4 

33.9-59.9 

Ph 

0 

,D 

a 

0.12 

0.148 

0.146 

0.119-0.176 

0.0954-0.212 

0.0808-0.242 


TO 

0.6 

0.776 

0.778 

0.638-0.922 

0.527-0.987 

0.461-0.997 

O 

v/s 

0.08 

0.0604 

0.0595 

0.0407-0.0801 

0.0244-0.101 

0.0149-0.116 


k-n 

0.5 

0.682 

0.698 

0.373-0.978 

0.228-1.08 

0.206-1.09 


Norm. 

9.9 

10.8 

10.9 

8.15-13.6 

6.40-14.8 

5.82-15.0 

£ 

A 

0.14 

0.199 

0.198 

0.132-0.267 

0.105-0.295 

0.101-0.299 

h-5 

TO 

0.6 

0.620 

0.602 

0.415-0.846 

0.302-0.975 

0.265-0.995 

v/ s 

0.20 

0.163 

0.162 

0.135-0.190 

0.121-0.208 

0.116-0.215 


k-n 

0.5 

0.651 

0.653 

0.347-0.955 

0.223-1.07 

0.205-1.09 


ables ((iV c h), ^ 2 ( 2 }, 773 ( 2 }). The V 3 samples have a somewhat 
larger variance than the others, in part due to the underlying 
noise in the model calculations and also because V 3 is explic¬ 
itly given a lower weight (recall that (iV c h) : 772 ( 2 } : ^2 {3} are 
weighted 1.2 : 1.0 : 0.6). 

The KLN results in the bottom row tell a somewhat dif¬ 
ferent story, as they cannot fit all observables simultaneously. 
While the fit to (N c h) is excellent, the ratio of 772 to V 3 is sim¬ 
ply too large and the model has no choice but to compromise 
between the two, similar to previous KLN results [39]. The 
posterior biases more towards V 2 than V 3 due to the explicit 
higher weight on V 2 . 

Figure [Tl] shows an expanded view of the 77 /s marginal dis¬ 
tributions for Glauber and KLN. The Glauber distribution is 
approximately normal with mean ~0.06 and 95% confidence 
interval ~ 0 . 02 - 0 . 10 , consistent with but mostly below 0.08. 
This is unsurprising and easily within the uncertainty of ex¬ 
isting results. KLN has a wider plateau-like distribution with 
mean ~0.16 and 95% confidence interval ~0.12-0.21. While 
the common estimate 0.20 was derived primarily from com¬ 
parisons to V 2 , the additional constraint from V 3 shifts the dis¬ 
tribution to somewhat smaller values and causes the plateau 
shape: rather than a strong peak, there is a range of values 
which all fit the data roughly equally. 

Table |U] quantitatively summarizes the posterior distribu¬ 
tions for each parameter including basic statistics, confidence 
intervals, and comparisons to previous estimates from earlier 
work with the same models EE1E1I1E1- All previous esti¬ 
mates fall within 95% confidence intervals, and most within 
68 %. 


V. CONCLUSION 

We have applied modern Bayesian methodology to system¬ 
atically compare an event-by-event heavy-ion collision model 
to experimental data. We chose a set of salient model param¬ 
eters including the shear viscosity 77 /s, evaluated the model 
over wide ranges of each parameter, and interpolated the re¬ 
sults with a Gaussian process emulator. Then, we used the 
emulator to calibrate the model to optimally reproduce exper¬ 
imental data and thereby extracted probability distributions 


for the true values of all model parameters simultaneously, 
including all correlations. 

When properly calibrated, the Monte Carlo Glauber model 
provides a good simultaneous fit to experimental multiplicity 
and flow data, while the Monte Carlo KLN model fails to 
simultaneously fit elliptic and triangular flow. The 77 /s distri¬ 
butions for the Glauber and KLN models are consistent with 
the commonly quoted values 0.08 and 0 . 20 , respectively, and 
in general the calibrated distributions reinforce and expand 
upon existing knowledge of these models. 

This study represents a significant step forward in state- 
of-the-art model-to-data comparison and establishes a frame¬ 
work for future analysis. Since the method does not reduce 
each parameter to a “best-fit” value but instead furnishes full 
probability distributions, it may be used to rigorously quan¬ 
tify uncertainties, examine correlations among parameters, 
and evaluate the efficacy of physical models, among other 
possibilities. It is easily extensible to arbitrary numbers of 
parameters and physical observables and to different models. 

Indeed, we plan to apply the methodology to a variety 
of other models, including the new initial condition model 
TrENTo— a flexible effective model which is ideal for this 
type of analysis [33]— and a 3+1D viscous hydrodynamics 
model with finite baryon chemical potential combined with 
recent data from the RHIC beam energy scan. By consid¬ 
ering data from multiple beam energies, we can probe the 
temperature dependence of 77 / s. 

We will include additional physical properties such as the 
size and shape of nucleons in the initial state, the hydro- 
dynamic equation of state, and the switching temperature 
from hydrodynamics to microscopic transport; and compare 
to more observables, e.g. identified particle spectra and dif¬ 
ferential flow. 

Finally, we anticipate upgrades to the methodology itself, 
notably more rigorous treatment of uncertainties and quan¬ 
tification of input-output correlations (analysis of variance). 
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APPENDIX: TRAINING THE EMULATOR 

Gaussian process emulators are non-parametric models 
(they do not assume a functional form) but they do require 
an assumed covariance function. One typically chooses a pa¬ 
rameterized functional form based on physical considerations, 
for example the squared-exponential function 

cr(x,x) = exp^- j 

generates smoothly-varying processes and is therefore com¬ 
patible with many models. If the model is known to oscillate, 
one would choose a periodic covariance function. Most mod¬ 
els have statistical noise which one accounts for by adding a 
diagonal noise term to the covariance function: 

a(x, x') ->• cr(x, x') + a^(5 xx /, 

where is the variance of the noise and <5 XX / is a Kronecker 
delta. 

Covariance functions often have free parameters known 
as hyperparameters , e.g. the squared-exponential correlation 
length £, which are not known a priori and must be estimated 
from model data. The selection of hyperparameters is known 
as training and may be accomplished by maximizing the like¬ 
lihood function H3 

logP(y|X,0) = -^y T E“ 1 y- i log |E y | - y log27r, (Al) 

where y is the vector of training outputs, X is the matrix 
of training input points, 6 is the vector of hyperparameters, 
and Ti y is the covariance matrix from applying the covariance 
function to the training data. Hence, the likelihood is the 
probability of observing the data given the model. The first 
term in the likelihood is the fit to data, the second term is 
a complexity penalty, and the third term is a normalization 
constant. 

This is best demonstrated by an example. The one¬ 
dimensional training data in Fig. [12] were generated from a 
Gaussian process with covariance function 

cr(x,x')= exp (- - ^ I ) + alS xx / (A2) 


Overfit: £ = 0.02, a n =0.001 



2 Oversmooth: £ = 3, a n =0.3 



2 Max. likelihood: £ = 0.462, a n =0.211 



-2 

0.0 


Actual: £ = 0.5, cr n =0.2 
T2 0 A T6 T8 


1.0 


x 


FIG. 12. Effect of varying the hyperparameters. Each panel 
shows a Gaussian process conditioned o n fab ricated training 
data using the covariance function Eq. (A2), where the line 
is the mean and the band is a 2 a confidence interval. The 
covariance function hyperparameters are different in each plot 
as indicated by the annotations. The data points are identical 
in each plot and were generated from a Gaussian process with 
hyperparameters annotated at the bottom. 


and hyperparameters 6 = (£, a n ) = (0.5, 0.2); now let us pre¬ 
tend we don’t know the true values of 6 and attempt to train 
a Gaussian process on the data. In the top panel of the figure, 
we use a short length scale £ with small noise cr n , so the GP 
interpolates each point exactly; however it rapidly wiggles and 
is almost certainly “overfit”. This choice of hyperparameters 
has a large complexity penalty and therefore a low likelihood. 
In the other extreme, we use a long length scale with large 
noise (middle panel), leading to a nearly linear GP that at¬ 
tributes most of the variance to noise. Here the likelihood 
is also low due to the poor fit to data. The most likely sce¬ 
nario is the compromise in the bottom panel, in which we 
estimate the hyperparameters by numerically maximizing the 
likelihood. Now, the trained GP smoothly interpolates the 
curvature of the training data while leaving some of the vari¬ 
ance as noise, true to the actual GP. 









14 


For the present study we use the covariance function given 
in Eq. (23) and restated here: 


/ /\ 2 

cr(x, x ) = cr GP exp 


-E 

L k= 1 


(xk - Xkf 

2 e k 


+ °n^x 


The hyperparameter ctq P is the overall variance of the Gaus¬ 
sian process and the £k are the independent length scales for 
each design parameter. We estimate the hyp erparameters by 
numerically maximizing the likelihood ( |A1| ) using a nonlin¬ 
ear conjugate gradient algorithm. Since the likelihood may 
have non-optimal local extrema, we repeat the optimization 
algorithm many times (minimum 100) from different random 
starting points. 

Table [In] lists the maximum-likelihood estimates of the hy¬ 
perparameters for each principal component in standardized 
units—input parameters scaled to [0,1] and principal compo¬ 
nents scaled to unit variance. We constrain the length scales 
to [0.3,10] for numerical robustness. 

In this work we fix the hyperparameters to the maximum- 
likelihood estimates during calibration. This neglects un¬ 
certainty in the hyperparameters themselves, although the 
present event-by-event model is well-behaved and the sample 
size is large, so varying the hyperparameters weakly affects 
the actual emulator predictions. But ideally, one would con¬ 
sider all predictions consistent with the data—not only the 
most likely—by MCMC-sampling the hyperparameter poste¬ 
riors. This significantly increases computational cost, since 
the GPs must be reconditioned for every set of hyperparam¬ 


eters, and conditioning requires computation of the inverse 
covariance matrix, an G(n 3 ) operation. We forgo this refine¬ 
ment until a future study. 


TABLE III. Maximum-likelihood estimates of the covariance 
function hyperparameters. 





Principal component 


1 

2 

3 

4 

5 


ctgp 

4.29 

2.20 

2.67 

0.923 

0.558 


i Norm 

1.89 

1.62 

1.95 

0.622 

0.300 

0 

£ a 

3.26 

1.61 

1.49 

0.579 

10.0 

1 

£ T 0 

1.20 

1.05 

1.65 

1.34 

0.300 

o 

£ rj/s 

1.69 

1.01 

1.17 

2.09 

10.0 


£ k-jr 

10.0 

4.77 

4.46 

10.0 

1.48 


O'n 

0.0349 

0.106 

0.558 

0.800 

0.933 


ctgp 

5.11 

3.36 

1.48 

1.28 

0.996 


£ Norm 

1.82 

1.39 

1.16 

0.907 

0.713 

£ 

£ A 

8.47 

4.54 

0.985 

1.10 

0.300 

J 

£ T 0 

0.927 

0.678 

0.808 

0.534 

0.359 


£ rj/s 

1.63 

0.851 

0.500 

0.434 

0.369 


£ ki r 

10.0 

8.33 

2.04 

1.43 

0.389 


T n 

0.0192 

0.0568 

0.807 

0.803 

0.606 
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